Plot Data
library(ggplot2)
# raw data
ggplot(dataset, aes(x=MMS, y=Counts)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_smooth(method=lm, formula = y ~ poly(x, 2), se=FALSE, colour="#666666") +
geom_point(aes(colour=Treatment, shape=Experiment), size=2) +
facet_grid(. ~ genotype) +
xlab(label = "MMS (%)") +
scale_shape_manual(values=15:19) +
scale_color_manual(values=c("#999999","#0072B2","#CC79A7","#009E73"))

# NormCounts Linear
ggplot(dataset[dataset$genotype %in% levels(dataset$genotype),], aes(x=MMS, y=NormCounts)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(colour="#333333") +
geom_smooth(method=lm, formula = y ~ x, se=TRUE, colour="#666666") +
facet_grid(. ~ genotype) +
xlab(label = "MMS (%)")

# NormCounts Quadratic
ggplot(dataset[dataset$genotype %in% levels(dataset$genotype),], aes(x=MMS, y=NormCounts)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(colour="#333333") +
geom_smooth(method=lm, formula = y ~ poly(x, 2), se=TRUE, colour="#666666") +
facet_grid(. ~ genotype) +
xlab(label = "MMS (%)")

# NormCounts Cubic
ggplot(dataset[dataset$genotype %in% levels(dataset$genotype),], aes(x=MMS, y=NormCounts)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(colour="#333333") +
geom_smooth(method=lm, formula = y ~ poly(x, 3), se=TRUE, colour="#666666") +
facet_grid(. ~ genotype) +
xlab(label = "MMS (%)")

# NormCounts2 Linear
ggplot(dataset[dataset$genotype %in% levels(dataset$genotype),], aes(x=MMS, y=NormCounts2)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(colour="#333333") +
geom_smooth(method=lm, formula = y ~ x, se=TRUE, colour="#666666") +
facet_grid(. ~ genotype) +
xlab(label = "MMS (%)")

# NormCounts2 Quadratic
ggplot(dataset[dataset$genotype %in% levels(dataset$genotype),], aes(x=MMS, y=NormCounts2)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(colour="#333333") +
geom_smooth(method=lm, formula = y ~ poly(x, 2), se=TRUE, colour="#666666") +
facet_grid(. ~ genotype) +
xlab(label = "MMS (%)")

# NormCounts2 Cubic
ggplot(dataset[dataset$genotype %in% levels(dataset$genotype),], aes(x=MMS, y=NormCounts2)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(colour="#333333") +
geom_smooth(method=lm, formula = y ~ poly(x, 3), se=TRUE, colour="#666666") +
facet_grid(. ~ genotype) +
xlab(label = "MMS (%)")

library(Cairo)
cairo_pdf("Figure3A_v1.pdf", width = 10, height = 4, family = "Arial")
datasubset <- dataset[dataset$genotype %in% levels(dataset$genotype),]
datasubset$genotype <- relevel(datasubset$genotype, ref = "WT")
ggplot(datasubset, aes(x=MMS, y=NormCounts2)) +
theme_bw() +
theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
axis.line = element_line(colour = "black"), text = element_text(size=14),
panel.border = element_blank(), panel.background = element_blank(),
axis.text.x = element_text(angle = 90, hjust = 1)) +
geom_point(aes(colour = genotype)) +
geom_smooth(method=lm, formula = y ~ x, se=TRUE, aes(colour = genotype), fill='#DDDDDD', size=0.5) +
facet_grid(. ~ genotype) +
xlab(label = "MMS (%)") +
ylab(label = "Normalized Counts") +
scale_x_continuous(labels = function(x) format(x, scientific = TRUE)) +
scale_color_manual(values=c('#000000','#000080','#808080','#800000'))
dev.off()
## quartz_off_screen
## 2
cairo_pdf("Figure3A_v2.pdf", width = 5, height = 4, family = "Arial")
datasubset <- dataset[dataset$genotype %in% levels(dataset$genotype),]
datasubset$genotype <- relevel(datasubset$genotype, ref = "WT")
ggplot(datasubset, aes(x=MMS, y=NormCounts2)) +
theme_bw() +
theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
axis.line = element_line(colour = "black"), text = element_text(size=14),
panel.border = element_blank(), panel.background = element_blank(),
axis.text.x = element_text(angle = 90, hjust = 1)) +
geom_point(aes(colour = genotype)) +
geom_smooth(method=lm, formula = y ~ x, se=TRUE, aes(colour = genotype), fill='#DDDDDD', size=0.5) +
#facet_grid(. ~ genotype) +
xlab(label = "MMS (%)") +
ylab(label = "Normalized Counts") +
scale_x_continuous(labels = function(x) format(x, scientific = TRUE)) +
scale_color_manual(values=c('#000000','#000080','#808080','#800000'))
dev.off()
## quartz_off_screen
## 2
Models
library(MASS)
library(DHARMa)
library(lme4)
library(lmerTest)
library(bbmle)
Linear formula
fit1 <- lm(Counts ~ Experiment + MMS*genotype, data = dataset)
print(summary(fit1))
##
## Call:
## lm(formula = Counts ~ Experiment + MMS * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -358.11 -116.30 48.29 106.75 364.56
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1925.11 100.06 19.239 < 2e-16 ***
## Experimentexp2 -343.00 69.07 -4.966 1.48e-05 ***
## Experimentexp3 -343.19 69.07 -4.969 1.46e-05 ***
## MMS -554124.81 152429.81 -3.635 0.000819 ***
## genotypeALC1 KO 319.15 129.79 2.459 0.018599 *
## genotypeALC1 KO PARP1 KO 134.34 129.79 1.035 0.307184
## genotypePARP1 KO -660.18 129.79 -5.087 1.01e-05 ***
## MMS:genotypeALC1 KO -841902.59 215568.30 -3.906 0.000373 ***
## MMS:genotypeALC1 KO PARP1 KO -687732.12 215568.30 -3.190 0.002848 **
## MMS:genotypePARP1 KO -108036.53 215568.30 -0.501 0.619142
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 195.4 on 38 degrees of freedom
## Multiple R-squared: 0.8914, Adjusted R-squared: 0.8657
## F-statistic: 34.66 on 9 and 38 DF, p-value: 1.242e-15
cat("AIC: ", AIC(fit1))
## AIC: 653.3867
simres <- simulateResiduals(fittedModel = fit1)
plot(simres)

fit2 <- lm(NormCounts ~ MMS*genotype, data = dataset)
print(summary(fit2))
##
## Call:
## lm(formula = NormCounts ~ MMS * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.10063 -0.02897 -0.00040 0.03175 0.08848
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.17580 0.02094 56.162 < 2e-16 ***
## MMS -370.11386 34.77268 -10.644 3.10e-13 ***
## genotypeALC1 KO 0.31434 0.02961 10.617 3.34e-13 ***
## genotypeALC1 KO PARP1 KO 0.29819 0.02961 10.071 1.57e-12 ***
## genotypePARP1 KO 0.25990 0.02961 8.778 7.20e-11 ***
## MMS:genotypeALC1 KO -661.77215 49.17600 -13.457 < 2e-16 ***
## MMS:genotypeALC1 KO PARP1 KO -627.76829 49.17600 -12.766 1.09e-15 ***
## MMS:genotypePARP1 KO -547.15757 49.17600 -11.127 8.14e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04456 on 40 degrees of freedom
## Multiple R-squared: 0.9843, Adjusted R-squared: 0.9816
## F-statistic: 359 on 7 and 40 DF, p-value: < 2.2e-16
cat("AIC: ", AIC(fit2))
## AIC: -153.1714
simres <- simulateResiduals(fittedModel = fit2)
plot(simres)

fit3 <- lm(NormCounts2 ~ MMS*genotype, data = dataset)
print(summary(fit3))
##
## Call:
## lm(formula = NormCounts2 ~ MMS * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.067414 -0.022126 -0.000339 0.021866 0.059274
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.985e-01 1.496e-02 66.727 < 2e-16 ***
## MMS -3.143e+02 2.485e+01 -12.646 1.47e-15 ***
## genotypeALC1 KO -1.661e-04 2.116e-02 -0.008 0.994
## genotypeALC1 KO PARP1 KO 8.597e-03 2.116e-02 0.406 0.687
## genotypePARP1 KO -9.859e-03 2.116e-02 -0.466 0.644
## MMS:genotypeALC1 KO -3.770e+02 3.515e+01 -10.726 2.46e-13 ***
## MMS:genotypeALC1 KO PARP1 KO -3.675e+02 3.515e+01 -10.455 5.26e-13 ***
## MMS:genotypePARP1 KO -3.173e+02 3.515e+01 -9.028 3.38e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03185 on 40 degrees of freedom
## Multiple R-squared: 0.9848, Adjusted R-squared: 0.9821
## F-statistic: 369.5 on 7 and 40 DF, p-value: < 2.2e-16
cat("AIC: ", AIC(fit3))
## AIC: -185.4122
simres <- simulateResiduals(fittedModel = fit3)
plot(simres)

fit4 <- lmer(Counts ~ MMS*genotype + (1|UID), data = dataset)
print(summary(fit4))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ MMS * genotype + (1 | UID)
## Data: dataset
##
## REML criterion at convergence: 470.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.65153 -0.51545 0.00911 0.43069 2.32751
##
## Random effects:
## Groups Name Variance Std.Dev.
## UID (Intercept) 69881 264.4
## Residual 11740 108.4
## Number of obs: 48, groups: UID, 12
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.696e+03 1.609e+02 9.088e+00 10.544 2.13e-06
## MMS -5.541e+05 8.455e+04 3.200e+01 -6.554 2.20e-07
## genotypeALC1 KO 3.192e+02 2.275e+02 9.088e+00 1.403 0.1939
## genotypeALC1 KO PARP1 KO 1.343e+02 2.275e+02 9.088e+00 0.590 0.5693
## genotypePARP1 KO -6.602e+02 2.275e+02 9.088e+00 -2.902 0.0174
## MMS:genotypeALC1 KO -8.419e+05 1.196e+05 3.200e+01 -7.041 5.54e-08
## MMS:genotypeALC1 KO PARP1 KO -6.877e+05 1.196e+05 3.200e+01 -5.752 2.23e-06
## MMS:genotypePARP1 KO -1.080e+05 1.196e+05 3.200e+01 -0.904 0.3730
##
## (Intercept) ***
## MMS ***
## genotypeALC1 KO
## genotypeALC1 KO PARP1 KO
## genotypePARP1 KO *
## MMS:genotypeALC1 KO ***
## MMS:genotypeALC1 KO PARP1 KO ***
## MMS:genotypePARP1 KO
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) MMS gALC1K gALKPK gPARPK MMS:AK MMSKPK
## MMS -0.250
## gntypALC1KO -0.707 0.177
## gALC1KOPARK -0.707 0.177 0.500
## gntyPARP1KO -0.707 0.177 0.500 0.500
## MMS:gALC1KO 0.177 -0.707 -0.250 -0.125 -0.125
## MMS:ALC1KPK 0.177 -0.707 -0.125 -0.250 -0.125 0.500
## MMS:PARP1KO 0.177 -0.707 -0.125 -0.125 -0.250 0.500 0.500
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
cat("AIC: ", AIC(fit4))
## AIC: 490.699
simres <- simulateResiduals(fittedModel = fit4)
plot(simres)

Quadratic formula
fit5 <- lm(Counts ~ Experiment + poly(MMS, 2)*genotype, data = dataset)
print(summary(fit5))
##
## Call:
## lm(formula = Counts ~ Experiment + poly(MMS, 2) * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -368.03 -127.99 45.33 111.29 361.68
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1661.90 72.86 22.810 < 2e-16 ***
## Experimentexp2 -343.00 72.86 -4.708 4.10e-05 ***
## Experimentexp3 -343.19 72.86 -4.710 4.07e-05 ***
## poly(MMS, 2)1 -1420.33 412.16 -3.446 0.001531 **
## poly(MMS, 2)2 65.59 412.16 0.159 0.874509
## genotypeALC1 KO -80.75 84.13 -0.960 0.343926
## genotypeALC1 KO PARP1 KO -192.33 84.13 -2.286 0.028601 *
## genotypePARP1 KO -711.50 84.13 -8.457 7.08e-10 ***
## poly(MMS, 2)1:genotypeALC1 KO -2157.97 582.88 -3.702 0.000753 ***
## poly(MMS, 2)2:genotypeALC1 KO 67.12 582.88 0.115 0.909001
## poly(MMS, 2)1:genotypeALC1 KO PARP1 KO -1762.80 582.88 -3.024 0.004719 **
## poly(MMS, 2)2:genotypeALC1 KO PARP1 KO -46.60 582.88 -0.080 0.936745
## poly(MMS, 2)1:genotypePARP1 KO -276.92 582.88 -0.475 0.637760
## poly(MMS, 2)2:genotypePARP1 KO -12.05 582.88 -0.021 0.983622
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 206.1 on 34 degrees of freedom
## Multiple R-squared: 0.8919, Adjusted R-squared: 0.8505
## F-statistic: 21.57 on 13 and 34 DF, p-value: 1.094e-12
cat("AIC: ", AIC(fit5))
## AIC: 661.1783
simres <- simulateResiduals(fittedModel = fit5)
plot(simres)

fit6 <- lm(NormCounts ~ poly(MMS, 2)*genotype, data = dataset)
print(summary(fit6))
##
## Call:
## lm(formula = NormCounts ~ poly(MMS, 2) * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.084112 -0.027962 -0.003658 0.033302 0.100678
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.000e+00 1.320e-02 75.760 < 2e-16
## poly(MMS, 2)1 -9.487e-01 9.145e-02 -10.374 2.31e-12
## poly(MMS, 2)2 3.251e-02 9.145e-02 0.355 0.724
## genotypeALC1 KO 2.243e-18 1.867e-02 0.000 1.000
## genotypeALC1 KO PARP1 KO 1.652e-17 1.867e-02 0.000 1.000
## genotypePARP1 KO 3.606e-17 1.867e-02 0.000 1.000
## poly(MMS, 2)1:genotypeALC1 KO -1.696e+00 1.293e-01 -13.116 2.80e-15
## poly(MMS, 2)2:genotypeALC1 KO 6.761e-02 1.293e-01 0.523 0.604
## poly(MMS, 2)1:genotypeALC1 KO PARP1 KO -1.609e+00 1.293e-01 -12.442 1.34e-14
## poly(MMS, 2)2:genotypeALC1 KO PARP1 KO -3.297e-02 1.293e-01 -0.255 0.800
## poly(MMS, 2)1:genotypePARP1 KO -1.402e+00 1.293e-01 -10.844 6.82e-13
## poly(MMS, 2)2:genotypePARP1 KO 4.241e-02 1.293e-01 0.328 0.745
##
## (Intercept) ***
## poly(MMS, 2)1 ***
## poly(MMS, 2)2
## genotypeALC1 KO
## genotypeALC1 KO PARP1 KO
## genotypePARP1 KO
## poly(MMS, 2)1:genotypeALC1 KO ***
## poly(MMS, 2)2:genotypeALC1 KO
## poly(MMS, 2)1:genotypeALC1 KO PARP1 KO ***
## poly(MMS, 2)2:genotypeALC1 KO PARP1 KO
## poly(MMS, 2)1:genotypePARP1 KO ***
## poly(MMS, 2)2:genotypePARP1 KO
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04572 on 36 degrees of freedom
## Multiple R-squared: 0.9852, Adjusted R-squared: 0.9806
## F-statistic: 217.2 on 11 and 36 DF, p-value: < 2.2e-16
cat("AIC: ", AIC(fit6))
## AIC: -147.7615
simres <- simulateResiduals(fittedModel = fit6)
plot(simres)

fit7 <- lm(NormCounts2 ~ poly(MMS, 2)*genotype, data = dataset)
print(summary(fit7))
##
## Call:
## lm(formula = NormCounts2 ~ poly(MMS, 2) * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.056350 -0.020282 -0.002504 0.023290 0.067449
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.849185 0.009453 89.834 < 2e-16
## poly(MMS, 2)1 -0.805602 0.065491 -12.301 1.87e-14
## poly(MMS, 2)2 0.027603 0.065491 0.421 0.676
## genotypeALC1 KO -0.179244 0.013368 -13.408 1.44e-15
## genotypeALC1 KO PARP1 KO -0.165958 0.013368 -12.414 1.43e-14
## genotypePARP1 KO -0.160591 0.013368 -12.013 3.73e-14
## poly(MMS, 2)1:genotypeALC1 KO -0.966345 0.092618 -10.434 1.97e-12
## poly(MMS, 2)2:genotypeALC1 KO 0.039466 0.092618 0.426 0.673
## poly(MMS, 2)1:genotypeALC1 KO PARP1 KO -0.941937 0.092618 -10.170 3.95e-12
## poly(MMS, 2)2:genotypeALC1 KO PARP1 KO -0.027920 0.092618 -0.301 0.765
## poly(MMS, 2)1:genotypePARP1 KO -0.813386 0.092618 -8.782 1.77e-10
## poly(MMS, 2)2:genotypePARP1 KO 0.023983 0.092618 0.259 0.797
##
## (Intercept) ***
## poly(MMS, 2)1 ***
## poly(MMS, 2)2
## genotypeALC1 KO ***
## genotypeALC1 KO PARP1 KO ***
## genotypePARP1 KO ***
## poly(MMS, 2)1:genotypeALC1 KO ***
## poly(MMS, 2)2:genotypeALC1 KO
## poly(MMS, 2)1:genotypeALC1 KO PARP1 KO ***
## poly(MMS, 2)2:genotypeALC1 KO PARP1 KO
## poly(MMS, 2)1:genotypePARP1 KO ***
## poly(MMS, 2)2:genotypePARP1 KO
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03275 on 36 degrees of freedom
## Multiple R-squared: 0.9855, Adjusted R-squared: 0.9811
## F-statistic: 222.6 on 11 and 36 DF, p-value: < 2.2e-16
cat("AIC: ", AIC(fit7))
## AIC: -179.8137
simres <- simulateResiduals(fittedModel = fit7)
plot(simres)

fit8 <- lmer(Counts ~ poly(MMS, 2)*genotype + (1|UID), data = dataset)
print(summary(fit8))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ poly(MMS, 2) * genotype + (1 | UID)
## Data: dataset
##
## REML criterion at convergence: 467.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.56492 -0.48396 -0.07458 0.48306 2.18606
##
## Random effects:
## Groups Name Variance Std.Dev.
## UID (Intercept) 69518 263.7
## Residual 13193 114.9
## Number of obs: 48, groups: UID, 12
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 1433.17 155.79 8.00 9.199
## poly(MMS, 2)1 -1420.33 229.72 28.00 -6.183
## poly(MMS, 2)2 65.59 229.72 28.00 0.285
## genotypeALC1 KO -80.75 220.33 8.00 -0.367
## genotypeALC1 KO PARP1 KO -192.33 220.33 8.00 -0.873
## genotypePARP1 KO -711.50 220.33 8.00 -3.229
## poly(MMS, 2)1:genotypeALC1 KO -2157.97 324.88 28.00 -6.642
## poly(MMS, 2)2:genotypeALC1 KO 67.12 324.88 28.00 0.207
## poly(MMS, 2)1:genotypeALC1 KO PARP1 KO -1762.80 324.88 28.00 -5.426
## poly(MMS, 2)2:genotypeALC1 KO PARP1 KO -46.60 324.88 28.00 -0.143
## poly(MMS, 2)1:genotypePARP1 KO -276.92 324.88 28.00 -0.852
## poly(MMS, 2)2:genotypePARP1 KO -12.05 324.88 28.00 -0.037
## Pr(>|t|)
## (Intercept) 1.58e-05 ***
## poly(MMS, 2)1 1.12e-06 ***
## poly(MMS, 2)2 0.7774
## genotypeALC1 KO 0.7235
## genotypeALC1 KO PARP1 KO 0.4081
## genotypePARP1 KO 0.0121 *
## poly(MMS, 2)1:genotypeALC1 KO 3.31e-07 ***
## poly(MMS, 2)2:genotypeALC1 KO 0.8378
## poly(MMS, 2)1:genotypeALC1 KO PARP1 KO 8.67e-06 ***
## poly(MMS, 2)2:genotypeALC1 KO PARP1 KO 0.8870
## poly(MMS, 2)1:genotypePARP1 KO 0.4012
## poly(MMS, 2)2:genotypePARP1 KO 0.9707
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p(MMS,2)1 p(MMS,2)2 gALC1K gALKPK gPARPK p(MMS,2)1:AK
## ply(MMS,2)1 0.000
## ply(MMS,2)2 0.000 0.000
## gntypALC1KO -0.707 0.000 0.000
## gALC1KOPARK -0.707 0.000 0.000 0.500
## gntyPARP1KO -0.707 0.000 0.000 0.500 0.500
## p(MMS,2)1:AK 0.000 -0.707 0.000 0.000 0.000 0.000
## p(MMS,2)2:AK 0.000 0.000 -0.707 0.000 0.000 0.000 0.000
## p(MMS,2)1KPK 0.000 -0.707 0.000 0.000 0.000 0.000 0.500
## p(MMS,2)2KPK 0.000 0.000 -0.707 0.000 0.000 0.000 0.000
## p(MMS,2)1:PK 0.000 -0.707 0.000 0.000 0.000 0.000 0.500
## p(MMS,2)2:PK 0.000 0.000 -0.707 0.000 0.000 0.000 0.000
## p(MMS,2)2:AK p(MMS,2)1KPK p(MMS,2)2KPK p(MMS,2)1:PK
## ply(MMS,2)1
## ply(MMS,2)2
## gntypALC1KO
## gALC1KOPARK
## gntyPARP1KO
## p(MMS,2)1:AK
## p(MMS,2)2:AK
## p(MMS,2)1KPK 0.000
## p(MMS,2)2KPK 0.500 0.000
## p(MMS,2)1:PK 0.000 0.500 0.000
## p(MMS,2)2:PK 0.500 0.000 0.500 0.000
cat("AIC: ", AIC(fit8))
## AIC: 495.3176
simres <- simulateResiduals(fittedModel = fit8)
plot(simres)

Cubic formula
fit9 <- lm(Counts ~ Experiment + poly(MMS, 3)*genotype, data = dataset)
print(summary(fit9))
##
## Call:
## lm(formula = Counts ~ Experiment + poly(MMS, 3) * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -363.06 -126.81 51.63 108.77 372.94
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1661.90 77.10 21.555 < 2e-16 ***
## Experimentexp2 -343.00 77.10 -4.449 0.000110 ***
## Experimentexp3 -343.19 77.10 -4.451 0.000109 ***
## poly(MMS, 3)1 -1420.33 436.14 -3.257 0.002798 **
## poly(MMS, 3)2 65.59 436.14 0.150 0.881472
## poly(MMS, 3)3 68.44 436.14 0.157 0.876361
## genotypeALC1 KO -80.75 89.03 -0.907 0.371616
## genotypeALC1 KO PARP1 KO -192.33 89.03 -2.160 0.038854 *
## genotypePARP1 KO -711.50 89.03 -7.992 6.4e-09 ***
## poly(MMS, 3)1:genotypeALC1 KO -2157.97 616.79 -3.499 0.001482 **
## poly(MMS, 3)2:genotypeALC1 KO 67.12 616.79 0.109 0.914069
## poly(MMS, 3)3:genotypeALC1 KO 128.29 616.79 0.208 0.836636
## poly(MMS, 3)1:genotypeALC1 KO PARP1 KO -1762.80 616.79 -2.858 0.007677 **
## poly(MMS, 3)2:genotypeALC1 KO PARP1 KO -46.60 616.79 -0.076 0.940275
## poly(MMS, 3)3:genotypeALC1 KO PARP1 KO 86.66 616.79 0.141 0.889202
## poly(MMS, 3)1:genotypePARP1 KO -276.92 616.79 -0.449 0.656682
## poly(MMS, 3)2:genotypePARP1 KO -12.05 616.79 -0.020 0.984537
## poly(MMS, 3)3:genotypePARP1 KO -110.34 616.79 -0.179 0.859227
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 218.1 on 30 degrees of freedom
## Multiple R-squared: 0.8932, Adjusted R-squared: 0.8326
## F-statistic: 14.75 on 17 and 30 DF, p-value: 3.574e-10
cat("AIC: ", AIC(fit9))
## AIC: 668.5997
simres <- simulateResiduals(fittedModel = fit9)
plot(simres)

fit10 <- lm(NormCounts ~ poly(MMS, 3)*genotype, data = dataset)
print(summary(fit10))
##
## Call:
## lm(formula = NormCounts ~ poly(MMS, 3) * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.069344 -0.027036 0.000241 0.025537 0.079450
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.000e+00 1.260e-02 79.381 < 2e-16
## poly(MMS, 3)1 -9.487e-01 8.728e-02 -10.870 2.84e-12
## poly(MMS, 3)2 3.251e-02 8.728e-02 0.372 0.712
## poly(MMS, 3)3 4.303e-02 8.728e-02 0.493 0.625
## genotypeALC1 KO 1.275e-16 1.782e-02 0.000 1.000
## genotypeALC1 KO PARP1 KO 1.334e-16 1.782e-02 0.000 1.000
## genotypePARP1 KO -8.574e-17 1.782e-02 0.000 1.000
## poly(MMS, 3)1:genotypeALC1 KO -1.696e+00 1.234e-01 -13.743 5.68e-15
## poly(MMS, 3)2:genotypeALC1 KO 6.761e-02 1.234e-01 0.548 0.588
## poly(MMS, 3)3:genotypeALC1 KO 1.309e-01 1.234e-01 1.061 0.297
## poly(MMS, 3)1:genotypeALC1 KO PARP1 KO -1.609e+00 1.234e-01 -13.037 2.40e-14
## poly(MMS, 3)2:genotypeALC1 KO PARP1 KO -3.297e-02 1.234e-01 -0.267 0.791
## poly(MMS, 3)3:genotypeALC1 KO PARP1 KO 9.865e-02 1.234e-01 0.799 0.430
## poly(MMS, 3)1:genotypePARP1 KO -1.402e+00 1.234e-01 -11.363 9.11e-13
## poly(MMS, 3)2:genotypePARP1 KO 4.241e-02 1.234e-01 0.344 0.733
## poly(MMS, 3)3:genotypePARP1 KO -1.146e-01 1.234e-01 -0.929 0.360
##
## (Intercept) ***
## poly(MMS, 3)1 ***
## poly(MMS, 3)2
## poly(MMS, 3)3
## genotypeALC1 KO
## genotypeALC1 KO PARP1 KO
## genotypePARP1 KO
## poly(MMS, 3)1:genotypeALC1 KO ***
## poly(MMS, 3)2:genotypeALC1 KO
## poly(MMS, 3)3:genotypeALC1 KO
## poly(MMS, 3)1:genotypeALC1 KO PARP1 KO ***
## poly(MMS, 3)2:genotypeALC1 KO PARP1 KO
## poly(MMS, 3)3:genotypeALC1 KO PARP1 KO
## poly(MMS, 3)1:genotypePARP1 KO ***
## poly(MMS, 3)2:genotypePARP1 KO
## poly(MMS, 3)3:genotypePARP1 KO
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04364 on 32 degrees of freedom
## Multiple R-squared: 0.988, Adjusted R-squared: 0.9823
## F-statistic: 175.4 on 15 and 32 DF, p-value: < 2.2e-16
cat("AIC: ", AIC(fit10))
## AIC: -149.8982
simres <- simulateResiduals(fittedModel = fit10)
plot(simres)

fit11 <- lm(NormCounts2 ~ poly(MMS, 3)*genotype, data = dataset)
print(summary(fit11))
##
## Call:
## lm(formula = NormCounts2 ~ poly(MMS, 3) * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.046457 -0.020582 0.000168 0.017199 0.054865
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.849185 0.009118 93.137 < 2e-16
## poly(MMS, 3)1 -0.805602 0.063168 -12.753 4.33e-14
## poly(MMS, 3)2 0.027603 0.063168 0.437 0.665
## poly(MMS, 3)3 0.036542 0.063168 0.578 0.567
## genotypeALC1 KO -0.179244 0.012894 -13.901 4.14e-15
## genotypeALC1 KO PARP1 KO -0.165958 0.012894 -12.871 3.39e-14
## genotypePARP1 KO -0.160591 0.012894 -12.455 8.18e-14
## poly(MMS, 3)1:genotypeALC1 KO -0.966345 0.089333 -10.817 3.20e-12
## poly(MMS, 3)2:genotypeALC1 KO 0.039466 0.089333 0.442 0.662
## poly(MMS, 3)3:genotypeALC1 KO 0.080002 0.089333 0.896 0.377
## poly(MMS, 3)1:genotypeALC1 KO PARP1 KO -0.941937 0.089333 -10.544 6.11e-12
## poly(MMS, 3)2:genotypeALC1 KO PARP1 KO -0.027920 0.089333 -0.313 0.757
## poly(MMS, 3)3:genotypeALC1 KO PARP1 KO 0.060256 0.089333 0.675 0.505
## poly(MMS, 3)1:genotypePARP1 KO -0.813386 0.089333 -9.105 2.13e-10
## poly(MMS, 3)2:genotypePARP1 KO 0.023983 0.089333 0.268 0.790
## poly(MMS, 3)3:genotypePARP1 KO -0.085854 0.089333 -0.961 0.344
##
## (Intercept) ***
## poly(MMS, 3)1 ***
## poly(MMS, 3)2
## poly(MMS, 3)3
## genotypeALC1 KO ***
## genotypeALC1 KO PARP1 KO ***
## genotypePARP1 KO ***
## poly(MMS, 3)1:genotypeALC1 KO ***
## poly(MMS, 3)2:genotypeALC1 KO
## poly(MMS, 3)3:genotypeALC1 KO
## poly(MMS, 3)1:genotypeALC1 KO PARP1 KO ***
## poly(MMS, 3)2:genotypeALC1 KO PARP1 KO
## poly(MMS, 3)3:genotypeALC1 KO PARP1 KO
## poly(MMS, 3)1:genotypePARP1 KO ***
## poly(MMS, 3)2:genotypePARP1 KO
## poly(MMS, 3)3:genotypePARP1 KO
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03158 on 32 degrees of freedom
## Multiple R-squared: 0.988, Adjusted R-squared: 0.9824
## F-statistic: 175.9 on 15 and 32 DF, p-value: < 2.2e-16
cat("AIC: ", AIC(fit11))
## AIC: -180.934
simres <- simulateResiduals(fittedModel = fit11)
plot(simres)

fit12 <- lmer(Counts ~ poly(MMS, 3)*genotype + (1|UID), data = dataset)
print(summary(fit12))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ poly(MMS, 3) * genotype + (1 | UID)
## Data: dataset
##
## REML criterion at convergence: 415
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.52884 -0.35459 -0.00389 0.40950 2.18082
##
## Random effects:
## Groups Name Variance Std.Dev.
## UID (Intercept) 69148 263.0
## Residual 14671 121.1
## Number of obs: 48, groups: UID, 12
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 1433.17 155.79 8.00 9.199
## poly(MMS, 3)1 -1420.33 242.25 24.00 -5.863
## poly(MMS, 3)2 65.59 242.25 24.00 0.271
## poly(MMS, 3)3 68.44 242.25 24.00 0.283
## genotypeALC1 KO -80.75 220.33 8.00 -0.367
## genotypeALC1 KO PARP1 KO -192.33 220.33 8.00 -0.873
## genotypePARP1 KO -711.50 220.33 8.00 -3.229
## poly(MMS, 3)1:genotypeALC1 KO -2157.97 342.59 24.00 -6.299
## poly(MMS, 3)2:genotypeALC1 KO 67.12 342.59 24.00 0.196
## poly(MMS, 3)3:genotypeALC1 KO 128.29 342.59 24.00 0.374
## poly(MMS, 3)1:genotypeALC1 KO PARP1 KO -1762.80 342.59 24.00 -5.145
## poly(MMS, 3)2:genotypeALC1 KO PARP1 KO -46.60 342.59 24.00 -0.136
## poly(MMS, 3)3:genotypeALC1 KO PARP1 KO 86.66 342.59 24.00 0.253
## poly(MMS, 3)1:genotypePARP1 KO -276.92 342.59 24.00 -0.808
## poly(MMS, 3)2:genotypePARP1 KO -12.05 342.59 24.00 -0.035
## poly(MMS, 3)3:genotypePARP1 KO -110.34 342.59 24.00 -0.322
## Pr(>|t|)
## (Intercept) 1.58e-05 ***
## poly(MMS, 3)1 4.78e-06 ***
## poly(MMS, 3)2 0.7889
## poly(MMS, 3)3 0.7800
## genotypeALC1 KO 0.7235
## genotypeALC1 KO PARP1 KO 0.4081
## genotypePARP1 KO 0.0121 *
## poly(MMS, 3)1:genotypeALC1 KO 1.64e-06 ***
## poly(MMS, 3)2:genotypeALC1 KO 0.8463
## poly(MMS, 3)3:genotypeALC1 KO 0.7113
## poly(MMS, 3)1:genotypeALC1 KO PARP1 KO 2.88e-05 ***
## poly(MMS, 3)2:genotypeALC1 KO PARP1 KO 0.8929
## poly(MMS, 3)3:genotypeALC1 KO PARP1 KO 0.8025
## poly(MMS, 3)1:genotypePARP1 KO 0.4269
## poly(MMS, 3)2:genotypePARP1 KO 0.9722
## poly(MMS, 3)3:genotypePARP1 KO 0.7502
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("AIC: ", AIC(fit12))
## AIC: 451.0196
simres <- simulateResiduals(fittedModel = fit12)
plot(simres)

Final Result
fit <- fit3
output <- coef(summary(fit))
output <- output[grep("MMS", rownames(output)),]
rownames(output) <- gsub("poly\\(|, [1-3]\\)","", rownames(output) )
rownames(output) <- gsub("genotype", paste0(" ",levels(dataset$genotype)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$genotype)[1], sep = " in " )
# suggested result table
kable(output, row.names = T)
| MMS in WT |
-314.2951 |
24.85330 |
-12.646013 |
0 |
| MMS: WT vs. ALC1 KO |
-377.0071 |
35.14787 |
-10.726313 |
0 |
| MMS: WT vs. ALC1 KO PARP1 KO |
-367.4845 |
35.14787 |
-10.455384 |
0 |
| MMS: WT vs. PARP1 KO |
-317.3321 |
35.14787 |
-9.028487 |
0 |
write.table(output, file = "Figure3A_Stats_Ref_WT.txt", quote = F, sep = "\t", row.names = T, col.names = NA)
# re-fit with ALC1KO reference
dataset$genotype <- relevel(dataset$genotype, ref = "ALC1 KO")
fit <- lm(NormCounts2 ~ MMS*genotype, data = dataset)
output <- coef(summary(fit))
output <- output[grep("MMS", rownames(output)),]
rownames(output) <- gsub("poly\\(|, [1-3]\\)","", rownames(output) )
rownames(output) <- gsub("genotype", paste0(" ",levels(dataset$genotype)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$genotype)[1], sep = " in " )
# suggested result table
kable(output, row.names = T)
| MMS in ALC1 KO |
-691.302188 |
24.85330 |
-27.8153108 |
0.0000000 |
| MMS: ALC1 KO vs. WT |
377.007067 |
35.14787 |
10.7263133 |
0.0000000 |
| MMS: ALC1 KO vs. ALC1 KO PARP1 KO |
9.522602 |
35.14787 |
0.2709297 |
0.7878383 |
| MMS: ALC1 KO vs. PARP1 KO |
59.674968 |
35.14787 |
1.6978260 |
0.0973087 |
write.table(output, file = "Figure3A_Stats_Ref_ALC1.txt", quote = F, sep = "\t", row.names = T, col.names = NA)